Fluctuation-driven Turing patterns 



Thomas Butler* and Nigel Goldenfeld 
Department of Physics and Institute for Genomic Biology, 
University of Illinois at Urbana Champaign, 1110 West Green Street, Urbana, IL 61801 USA 

(Dated: February 2, 2011) 

Models of diffusion driven pattern formation that rely on the Turing mechanism are utilized in 
many areas of science. However, many such models suffer from the defect of requiring fine tuning 
of parameters or an unrealistic separation of scales in the diffusivities of the constituents of the 
system in order to predict the formation of spatial patterns. In the context of a very generic 
model of ecological pattern formation, we show that the inclusion of intrinsic noise in Turing models 
leads to the formation of "quasi-patterns" that form in generic regions of parameter space and 
are experimentally distinguishable from standard Turing patterns. The existence of quasi-patterns 
removes the need for unphysical fine tuning or separation of scales in the application of Turing 
models to real systems. 
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The study of the emergent spatiotemporal patterns in 
physical or biological systems is an exciting and fruit- 
ful line of research in physics and in many other dis- 
ciplines such as chemistry, ecology, animal biology, and 
neuroscience [IHS] . Examples include patterns on animal 
coats [5J, engineered bacterial systems [7], chemical pat- 
tern formation [8j, mussel population densities [5], and 
Rayleigh-Benard convection in fluids [TO] . 

One particularly satisfying aspect of these studies is 
that insight into the origins of one kind of pattern often 
yields insight into the origins of patterns in entirely dif- 
ferent systems. A key example is the Turing mechanism 
[3J. Turing's argument, which will be described in detail 
below, showed how diffusion, which is typically thought 
of as a randomizing influence, can give rise to spatial 
pattern formation when there are two or more classes of 
degrees of freedom (species) with "activator" and "in- 
hibitor" dynamics. This mechanism has been proposed 
as an explanation for an enormous variety of systems in- 
cluding short (< 10m) length scale patchiness in plank- 
tonic ecosystems [TTlfH] . patterning in plant-resource 
systems [TS], patchiness in insect abundance [TB], stripe 
and spot patterns on the coats of animals [BJ , patterns in 
mussel beds [9] and even the geometric visual hallucina- 
tions experienced by shamans and users of hallucinogenic 

drags gHHI- 
However, in spite of the seeming success of the Tur- 
ing mechanism in explaining patterns across many disci- 
plines, the partial differential equations representing the 
dynamics of systems with Turing patterns typically re- 
quire unphysical fine tuning of parameters or separation 
of scales in the diffusivities of the different species in or- 
der to predict pattern formation [5[ l8[ HT | [T8H2T] . The 
requirement that the system either have fine tuning of 
kinetic parameters or a separation of scales in diffusivi- 
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ties in order to predict patterns, is unphysical for many 
applications. Both the of these issues will be referred 
to below as the "fine tuning problem," even though fine 
tuning is only strictly needed when a separation of scales 
in diffusivities is not present. To resolve the fine tuning 
problem for Turing patterns we show that a full statistical 
mechanical treatment of Turing patterns, where fluctua- 
tions due to the discrete nature of the degrees of freedom 
in the system - intrinsic noise - are included, the fine 
tuning problem is resolved [2"T] , 

It may seem counterintuitive to claim that including 
fluctuations resolves the fine tuning problem for Turing 
patterns because fluctuations are generally expected to 
destabilize ordered states such as spatial patterns. This is 
the rule in standard statistical mechanics 22! and many 
statistical mechanical models in ecology [33J HI] ■ How- 
ever, exceptions exist in systems out of equilibrium. For 
example, careful experiments on Rayleigh-Benard con- 
vection have shown that fluctuations can drive the for- 
mation of convection rolls in fluid dynamics that would 
not form in the absence of fluctuations [25]. In ecology, 
recent theoretical work and careful data analysis have 
shown that the observed cyclic population dynamics of 
predator-prey systems can be explained in many cases 
by fluctuation driven cycles in time |26H29j . Similar phe- 
nomena have been predicted in evolutionary game theory 
and systems biology [301 |3T] • In cell biology, simulating 
the interactions of individual proteins in discrete time 
and space in a model of proteins that regulate cell divi- 
sion in e-coli results in pattern formation over a wider 
range of parameters than the corresponding reaction- 
diffusion partial differential equations |32j . Thus it seems 
possible that a full many body treatment of the Turing 
mechanism that incorporates intrinsic noise will resolve 
the fine tuning problem. 

The purpose of this paper is to present an analyisis 
of the Turing mechanism with intrinsic noise included to 
resolve the fine tuning problem. The analysis results in a 
derivation of a phase diagram and to power spectra with 
experimentally distinctive and relevant properties. This 
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paper is an expansion and elaboration of our paper [5T] 
which originally reported the resolution of the fine tuning 
problem of Turing instabilities through the incorporation 
of intrinsic noise. We will first review the Turing mecha- 
nism, and then present an extremely simple model of the 
Turing mechanism for planktonic predator-prey popula- 
tions that we then analyze in detail. The results of the 
analysis show that in large regions of parameter space 
predicted by deterministic modeling to have only triv- 
ial spatial states a new kind of spatial pattern emerges 
that we call a "quasi-pattern." The quasi-pattern state is 
analogous to intrinsic noise driven "quasi-cycles" recently 
discovered in the time domain |26j . Quasi-patterns are 
recognizable immediately as spatial patterns, but with a 
few important, experimentally relevant, differences from 
patterns predicted with deterministic analysis. The fi- 
nal sections of the paper will focus on possible experi- 
mental tests and extensions of the theory developed in 
the body of the paper. We focus on a model of plank- 
tonic predator-prey interactions throughout the paper for 
simplicity and also because predator-prey systems have 
been extensively analyzed theoretically [m[T21[T51l2"tJll55] 
and there is beginning to be an experimental literature 
[T5] . However, we emphasize that the goal of this 
paper is insight into the general interactions of intrin- 
sic fluctuations with the Turing mechanism for pattern 
formation and that the results should be valid for most 
models of Turing instabilities. Evidence for this asser- 
tion is provided by the recent replication of our results on 
the Brusselator model of chemical pattern formation |34j , 
which additionally pointed out that the phase boundary 
can differ for different species of reactants, a model of 
embryonic pattern formation [35] . as well as our own 
forthcoming results on pattern formation on the visual 
cortex [17] , 



I. THE TURING MECHANISM 

The Turing mechanism in its most basic form requires 
two different species that react and diffuse. One species, 
the "activator," diffuses relatively slowly, and catalyzes 
(activates) both its own production and the production of 
the second species. The second species, the "inhibitor," 
diffuses faster, and reduces (inhibits) the concentration 
of both the activator species and itself. These combined 
mechanisms lead to pattern formation from random ini- 
tial conditions. We illustrate the mechanism with the 
example of predator-prey dynamics with random initial 
conditions 

1. Random regions of activator (prey) with higher lo- 
cal concentrations reproduce rapidly, leading to dense 
clumps of activator species that then begin to diffuse. 

2. Rapidly diffusing inhibitors (predators) are produced 
in the neighborhood of the high density autocatalyzing 
clumps of prey. 

3. The predators inhibit the spread of the prey clumps 
through their production in the neighborhood of prey 



clumps. The autoinhibitory nature of predators prevents 
them from overwhelming the prey population. 

These steps, summarized in fig. [T] show how activator- 
inhibitor dynamics can lead to spontaneous pattern for- 
mation [3 . As was noted above, formalizing this argu- 
ment into standard deterministic reaction-diffusion equa- 
tions results in models that only exhibit Turing patterns 
if the predator (inhibitor) diffusivity is much larger than 
the prey (activator) diffusivity or the parameters are fine 
tuned [3ll8llTTl[T8H2"0]. Note that consistent with the ex- 
istence of pattern forming systems which do not appar- 
ently display very large separation of diffusivities [T51 [TS] 
the qualitative argument made above for pattern forma- 
tion does not depend on very large differences in diffu- 
sivities, nor on additional kinetic details. 




FIG. 1: Illustration of the steps of the Turing mechanism as 
described in the text. The figure should be viewed from top to 
bottom. The prey (activators) are represented by black dots, 
and the predators (inhibitors) are represented by red dots. 



II. TURING PATTERNS IN THE LEVIN-SEGEL 
MODEL 

One of the simplest models of Turing patterns is drawn 
from ecological pattern formation and was originally 
introduced to model plankton- herbivore dynamics [llj . 
The reaction diffusion equations for this model are 

d t tp = ^V 2 -0 + bip + eij} 2 - (pi + P2)4"-P 
d t ip — v\I 2 Lp + P2<ptp — dip 2 
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where the plankton population ip is the activator, as can 
be seen by the nonlinear growth term eip 2 , and the her- 
bivore population ip is the inhibitor due to the predation 
terms pipip and the competition term —dip 2 . The non- 
linear growth term eip 2 was origingally introduced to be 
a proxy for predator satiation but can also be in- 
terpreted as an Allee effect, wherein many species have 
enhanced reproduction at higher concentrations (for a 
review, see [3l)j). 

Setting pi = and P2 — p, the model contains a stable 
homogeneous coexistence state when 



p > e and p 2 > de 
with stationary fixed point populations given by 



Ips = 



bd 



bp 



p 2 — de 7 p 2 — de 

It contains a Turing instability if 



(2) 



(3) 



- > 



fl \(\/p/d- Vp/ d - e /p) 



(4) 



The model is only valid when the coexistence fixed point 
is stable. Outside of that regime, a plankton regulation 
term such as —fip 3 is required to make the model valid. 
For the present analysis, we assume that fip 3 is suffi- 
ciently small to be ignored. The fine tuning problem can 
be illustrated in this model by taking the set of 0(1) ki- 
netic parameters 6=1/2, e — 1/2, d = 1/2 and p = 1. 
With these parameters Eq. [4] shows that non-generic dif- 
fusivitics, v/fi > 27.8, are required for pattern formation. 
Similar results are obtained for other generic parameter 
sets. 



A. Extrinsic noise driven pattern formation 

To gain preliminary insight into the effects of intrinsic 
noise on the Levin-Segel model, we can analyze the effects 
of unconserved extrinsic noise on the linearized dynamics 
of the Levin-Segel model as in previous studies of extrin- 
sic noise driven pattern formation [2TJ . While not 
identical, expansion schemes such as the system size ex- 
pansion [JO] indicate that the effects of extrinsic noise 
and intrinsic noise on the linearized dynamics of reaction 
diffusion systems are closely related. Additionally, we 
will use the calculation of the effects of extrinsic noise on 
the Levin-Segel model to predict observable differences 
between unconserved extrinsic and intrinsic noise driven 
pattern formation. 

To calculate the effects of extrinsic noise, we look at the 
Fourier transformed dynamics of the fluctuations from 
the coexistence fixed point with added white noise £, vari- 
ance C. These dynamics are given by 



The matrix A is the Fourier transformed stability ma- 
trix and x is the vector of deviations from equilibrium of 
predator and prey populations respectively, 



—vk 2 — pips 

-Pips 



Pfs 

—\xk 2 + eip s 



(6) 



Simple manipulations yield the average power spectrum 



P(k,u) = C [pVs + (e*Ps ~ M 2 ) 2 ] x 



(pbip s + fivk 4 — uj 2 



^ ev (i -Pt)f+ w 2 ((e - P )yj s - (m + v)k 2 ) 
\ ev / 



1 -2 



(7) 



iojx = Ax + £ 



(5) 



To a crude approximation, Eq. [7] predicts that pat- 
terns (indicated by peaks in the power spectrum) form 
whenever ev > p[i, and that without noise and away from 
a classical Turing instability the power spectrum is zero. 
As anticipated, the condition ev > pjj, can be satisfied 
easily and avoids the fine tuning problem. However, the 
calculation with the extrinsic noise considered here dif- 
fers in important ways from the intrinsic noise case, such 
as the determination of the strength of the noise and the 
presence of diffusive noise. As will be shown below, these 
differences lead to experimentally distinguishable differ- 
ences in the resulting spatiotemporal patterns. 



III. PREDATOR-PREY MODEL WITH 
INTRINSIC NOISE 

To systematically include the effects of intrinsic noise 
requires a model defined at the level of individual organ- 
isms, since intrinsic noise is generated by the stochas- 
tic nature of individual birth and death events as well 
as the stochastic interactions between individual organ- 
isms. Such a description of the dynamics at the individ- 
ual level is called an individual level model (ILM). One 
simple way to define an ILM is to specify the reactions 
that can take place in a well mixed patch of volume V. 
To include space, a lattice of patches can be considered 
with additional reactions corresponding to movement of 
predator and prey organisms between the patches. With 
parameters to specify the relative rates of the reactions, 
a model of individual level interactions on a single patch 
that incorporates intrinsic noise is fully specified. 

For an ILM version of the Levin-Segel model we con- 
sider the following reactions 

P\ PP 

e/V 

pp U ppp 

PH P1 4 V H 

phH v hh 

HH d ^H (8) 
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where P denotes plankton and H denotes herbivores, 
with the parameters as described above. Stochastic tra- 
jectories of H and P, enumerated by m and n respec- 
tively, are described by the master equation 

d t P(m, n) — b(—nP(m, n) + (n — l)P(m, n — 1)) 
+ ~((n- l)(n-2)P(m,n- 1) - n(n - l)P(m,n)) 

+ y(-mriF(m, n) + (m)(n + l)P(m, n + 1)) 

P2 

+ — (— mnP(m, n) + (m — l)(n + l)P(m — 1, n + 1)) 

+ — [(m + l)mP{m + 1, n) — m(m — l)P(m, ri)] (9) 

The master equation, which is exactly equivalent to 
the specification of the model as a collection of reactions 
in Eq. [H] can then be used to analyze the ILM version 
of the Levin-Segel model by applying techniques from 
non-equilibrium statistical mechanics. 

A. Field theory representation of the model 

While several options exist for analysis of the master 
equation, such as direct expansion of the master equa- 
tion |40j . we analyze the master equation by a mapping 
to field theory, because it is convenient for handling spa- 
tially extended systems. To analyze the master equation 
using the techniques of field theory, we introduce the op- 
erators 

a\m, n) = m\m — 1, n) 
a\m, n) — \m + 1, n) 

[a, a] = 1 
c|m, n) = n\m, n — 1) 
c|m, n) = \m, n + 1) 

[c,c] = l (10) 

and the state = ^2 P(n)\n). These definitions allow 
the master equation to be mapped to a bosonic field the- 



ory [4"TM45| . As an explicit example of how to convert 
the master equation to a field theory, consider the mas- 
ter equation corresponding to the second reaction in Eq. 
[8] alone. 

d t P{n) = y[(n-l)(n-2)P(n-l)-n(n-l)P(ri)] (11) 

Ignoring V for now, we multiply both sides by \n) and 
sum over n 

J2dtP(n)\n} = 

e \{n - l)(n - 2)P(n - 1) - n(n - l)P(n)] |n) (12) 

n 

We next shift the sums, and manipulate the first term 
in the sum 

Let n' = n — 1 — s- n = n' + 1 . 

e^V(n' - l)P(n')K + 1) ,n' -> n 

n' 

= ec 3 c 2 J2P(n)\n) 

= ec 3 c 2 |^) (13) 
We now work out the second term in the sum 
e n(n — l)P(n) \n) 

n 

= ecccc\ip} (14) 

This yields 

d t |V}=e[c 3 -c 2 ]c 2 |^> (15) 

Similar analyses lead to second quantized forms for 
the rest of the master equation. We can now assemble 
the entire Hamiltonian. We start by writing the master 
equation in second quantized form 



d t W) 



&(c 2 -c)c+^(c 3 



Pi 



c )c + — (aac 



— (a ac — aacc) 



— (1 -a)aa 2 



(16) 



Since the standard definition of the Hamiltonian is 

d t \i>) = -H\$) (17) 

we have 

—H = b(c 2 — c)c + y(c 3 ~ c 2 )c 2 + ^(aac — aacc) 
+ ^-{a 2 ac — aacc) + — (1 — a)aa 2 (18) 



According to the standard mapping using coherent 
states between Hamiltonians represented by bosonic op- 
erators and functional integral representations of the 
same dynamics with Lagrangians, we can write down 
the Lagrangian, generalized to space. As in quantum 
mechanics, the mapping can be worked out for general 
Hamiltonians [HI US]. To generalize to space, we im- 
plement a random walk between patches of volume V 
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for every organism as a reaction with rate tj where i 
is an index for species. Appropriately rescaled |28j . the 
continuum limit and mapping to the functional integral 
formulation yields the Lagrangian 

C = adta + cdtc — vd\7 2 a — /ic\7 2 c 

+H(c,a 7 c,a) (19) 



In the Lagrangian formulation of Eq. 19 a, b and 
their conjugate variables are no longer operators, but are 
functions that are integrated over as in standard bosonic 
functional integrals. The starred variables loosely corre- 
spond to noise and the unstarred to values of predator 
and prey, but direct physical interpretation is not trivial 
[4"Tl H5] . The initial conditions are ignored, because the 
focus of this paper is the long time limit and there is only 
one attractor in the system. 

To transform to more physical variables, the standard 
Cole-Hopf transformation can be applied to transform 
the field variables to direct number and noise represen- 
tations. This transformation is given by 



ze 



c = pe 9 



(20) 
(21) 



the new field variables z and p can be heuristically inter- 
preted as the number of predator and prey respectively 
(the precise interpretation is that their expectation val- 
ues correspond, i.e. (f(p,z)) = (f(Np,Nii))) and the 
auxiliary fields denoted by carets generate the intrinsic 
noise, as will be seen below by showing that the minimum 
of the action, which corresponds to mean field theory is 
at p = z — 0. The Lagrangian in the new variables is 

C = xdt.z + pdtp — vz\7 2 z — pp\7 2 p 

- uz{Vz) 2 - pp{Vp) 2 + bp{l - e p ) 
Pi 



+ - e") + ^zp{\ - e-") 



zp(l 



V 

3^) + ^ 2 (l-e- £ ) 



(22) 



B. System size expansion 



We now can carry out the system size expansion in 
the field theoretic formalism. Other than notation, it is 
identical to the direct expansion of the master equation 
reviewed in [10]. We expand the fields as 



V 



P 



Vtp + VVt], p = Vip + 



P_ 



(23) 



To perform this expansion to consistent order, it is 
necessary to expand the exponentials out to second order. 



This is because the expansion will promote second order 
terms to first order. The result is an expansion of the 
Lagrangian in the form 



C = Wd+ C 2 + 0(1/VV) 



(24) 



We once again carry out the expansion explicitly for 
the term coupled by e/V. 



(Vv> + Wz) (Vv> + Wz) 



i- i 



-< ! Wi; 2 p + + 2£p ) + 0(1/ W) 



VV 2V 
(25) 



Collecting terms of leading order, yV, we have 

C\ = pdtip + zdt(p — vzV 2 (p — /ipV 2 ?/; — bipp — eip 2 p 
+ btpipp — cipip(z — p) + dip 2 z (26) 

It is trivial to extract the mean field PDE's by using 
the Euler-Lagrange equations. The equations that result 
are 



~6z~ 



d t ip - pV 2 ip + bip + etp 2 - (pi + p 2 )ipip = 



(27) 



which is the first of the equations for the Levin-Segel 
model. The second equation is 



= d t tp - vV 2 ip + p 2 ipip -dip 2 = (28) 



again reproducing the Levin-Segel model equation of mo- 
tion. Note that the auxiliary fields have zero expectation 
value at mean field, which confirms the interpretation 
that they correspond to noise. Now £ 2 can be assembled. 
The terms in d that are linear in r\ or £ correspond to the 
stability matrix of the MFT. Terms that are quadratic in 
the hatted variables p and z are noise terms and will be 
considered next. 

Proceeding, we have 

d = zd t i] + pdti - zv\I 2 j] - ppSI 2 i + Pi'qipp 
— Pmi>{z — p) + 1di}ipz + b£p + 2e£'i/jp + b^ipp 
-P2&&-P) (29) 

We convert this into a Fourier transformed matrix form 
that includes time and space 



d = V T d t x - y T Ax - \) T By 



with vectors given by 



' \\)- V= \p 



(30) 



(31) 



so that the predator variables are on top. The matrix A 
is the Jacobian of the MFT with space and is given by 

A _ / -uk 2 +p 2 ip-2d<p p 2 <p \ fiO) 

~ V -(P1+P2W ~p.k 2 +b+2e,p~(p 1 +p2)<p / ^ ' 

The matrix for the correlations of the noise is given by 

B = ( d P 2 +P2'PTf' + l "Pk 2 -p2<Pi> \ /OOj 

\ —P2<pip btp+eip 2 +bip^} J rP2^pip+^Lipk 2 J ^ ' 

We also now note that C 2 is in the form of a Lagrangian 
in the Martin-Siggia-Rose (MSR) response function for- 
malism for Langevin equations [121 [50] . 



C. The power spectrum 

We now extract the stochastic differential equations 
(SDE) that govern the dynamics of the fluctuations, and 
calculate the power spectrum of the fluctuations. The 
Langevin equations from the response function formalism 
are 



iuix = Ax + 7(0;) 
(7i(w)7j(-^)) = B ij 



We solve formally to obtain 

x = (A + iuj)- 1 ^^) = D(uj)- 1 -y(u)) (35) 
The power spectrum is 

(p 22 7i - D 12l2 )(D* 22ll - D* 12l2 )) 



(XiX*) 



\det(D)\ 2 

\D 22 \ 2 B n - 2D 12 Re{D 22 )B 21 + \D 12 \ 2 B 2 



\det{D)\ 2 



(36) 



(34) 



To find the phase diagram, take p\ — > 0, p 2 = p. This 
simplification does not substantially change the dynamics 
of the model. In terms of elements of the stability matrix 
J = A(k = 0, u> = 0), the denominator of the power 
spectrum is 

det(D) = (Jn + iu - vk 2 )(J 22 + ito - ^k 2 ) - J\ 2 J 2 \ 
= det(J) + ioj(Tr(J) - (/i + v)k 2 ) 
- (Jn/i + J 22 v)k 2 + ^fc 4 - uj 2 (37) 

The full expression for the power spectrum is 



P(k,u) 



\D 22 \ 2 B n - 2D 12 Re(D 22 )B 21 + \D 12 \ 2 B 22 



(det{J) + [iisk 4 -uj 2 - {Juiik 2 + J 22 vk 2 )) 2 + w 2 (Tr(J) - (p + v)k 2 f 



(38) 



Recall the fixed point values at coexistence are Now we evaluate the determinant of the ODE stability 

matrix (J above) and the trace 

pb 
p A — de 

db det ( J ) = P# ( 41 ) 

^ = ~2 T 

p A — de 

(39) The trace is 
Using the fixed point values, the matrix A can be further 

simplified to Tr(J) = (e-p)ijj (42) 



* ^ —pip -{ik 2 + eip ^ ! '''' 



Simplifying the denominator in Eq|38| yields 



\det(D)\ 2 = (det(J) + ^vk 4 - uj 2 - [J llt ik 2 + J 22 vk 2 )) 2 + uj 2 (Tr(J) - (/z + v)k 2 ) 2 
= (pbip + \ivk 4 - uj 2 - tpi-p^ik 2 + evk 2 )) 2 + uj 2 ((e - p)tp - (fj, + v)k 2 ) 2 

= (pbip + [ivk 4 - lo 2 - ^k 2 ev (l - ^) ) 2 + u 2 ((e - p)ip - (fx + v)k 2 ) 2 (43) 



The form of the denominator for lu — is (A — Bk 2 + Ck ) 2 , which has a minimum at non zero k. This mini- 
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mum corresponds to an emergent length scale, and is the Notice the symmetry in the noise correlations. We now 
first indication of pattern formation. Systematic demon- can expand out the numerator of Eq. 38 
stration of the emergence of pattern formation requires 
accounting for the k dependence in the numerator. The 
noise matrix B can be simplified to 

B= { Zpf^P + vpk 2 -piptp \ _ 
I —pipip 2ptpip + /iipk 2 J ^ ' 



\D 22 \ 2 B n - 2D 12 Re{D 22 )B 21 + \D 12 \ 2 B 22 

= \eij) — /ifc 2 + iuj\ 2 (2ptptp + vipk 2 ) + 2pip(eip — /J,k 2 )(pp , ip) + p 2 (p 2 (2ptpip + /itpk 2 ) 

= (eip - \ik 2 ) 2 {2pipi\) + vpk 2 ) + oj 2 (2pipip + vpk 2 ) + 2p 2 <p 2 tp(eip - fik 2 ) + p 2 <p 2 (2ppip + ^tpk 2 ) (45) 



This gives the final form of the power spectrum 



P{k,u) 



(eip - /ik 2 ) 2 (2ptpip + vipk 2 ) + Lo 2 (2p(piJj + vipk 2 ) + 2p 2 ip 2 Tp(e4> - pik 2 ) + p 2 Lp 2 (2ptpip + fiipk 2 ) 



(pbip + (Mfk 4 - uj 2 - %pk 2 ev (l - + cj 2 ((e - p)ip -(fi+ v)k 2 ) 2 



(46) 



IV. ANALYSIS OF THE POWER SPECTRUM 



A. Phase diagram for quasi-patterns 



typical of generic parameters. The phase diagram of the 
system bears out this conclusion as shown in figure [2] 



The expression for the power spectrum in Eq. [46] is 
not very illuminating, and does not simplify a great deal. 
To find quasi-patterns we note that the highest power of 
k in the denominator of Eq. [46] is larger than the highest 
power in the numerator. That means for sufficiently large 
k, the power spectrum is always decreasing. Thus, to 
show the existence of a maximum, it is sufficient to show 
that for small k, the power spectrum is increasing. This 
can be shown by computing and evaluating at k 2 = 0. 
When this expression is greater than 0, there is pattern 
formation. This yields the analytical criterion 



> 



p 3 (5p 2 + 7de) 



PL e(4p 4 + 5p 2 de + M 2 e 2 



(47) 



This criterion is much less stringent than the criterion 
for Turing instabilities. The conditions for a Turing in- 
stability are 



'- > h ' 

M V y/p/d - ^p/d - e/p 



(48) 



For the generic parameters 6 = 1/2, p = 1, d = 
1/2, e = 1/2 the criterion [47] yields v/fj. > 2.48, while 
the Turing condition yields vf[i > 27.8. This behavior is 




FIG. 2: Phase diagram over stable parameter region in p/d. 
Region I is MFT level pattern formation, region II contains 
fluctuation driven quasi-patterns, and region III is a spatially 
homogeneous phase. 



An additional feature of the model is that oscillations 
and spatial pattern formation are essentially decoupled. 
This means that the model predicts global population 
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oscillations and spatial pattern formation, but not trav- 
eling waves. The mathematical origin of this can be seen 
in Eq. [7] The k 2 term with a negative coefficient at 
uj = is quickly overwhelmed by the positive k 2 depen- 
dence of the lo 2 term as the frequency begins to grow. In 
the power spectrum (fig. [3]) this can be seen as the deep 
valley between the peaks in k and oj. This interpretation 
is supported by preliminary simulations. 




k 



FIG. 3: Power spectrum with p=l, vj\i=\h 



B. Wavelength of fluctuation driven patterns 

To a fairly good approximation, the wavelength of the 
Turing quasi-patterns can be calculated. The wavelength 
corresponds to the wave vector that maximizes the power 
spectrum. To calculate that value, consider the numera- 
tor of the power spectrum only at u> = 0. 

(pbip + ^k 4 ~ ^k 2 ev (l - ^) ) V (49) 

The minimum of this expression will correspond with 
reasonable accuracy to the real wavelength and can be 
obtained through straightforward calculation to be 

X m = ^= (50) 

k m V W \ ev ' 

This shows that for a fixed ratio of diffusivities, the 
wavelength increases as the square root of the diffusivity. 
In addition, while the phase diagram of the system (fig 
|2j) and therefore the presence of Turing quasi-patterns 
depends on diffusivity only through the ratio v/n, the 
wavelength of the patterns depends on the values of the 
diffusivities. 

This calculation also implies that the wavelength of 
the quasi-patterns is closely related to the wavelength of 



patterns in the region of the phase diagram where pat- 
terns are generated at mean field. In the standard theory 
of Turing patterns, patterns are formed when the homo- 
geneous steady state is unstable to perturbations with a 
specific set of wave vectors k. The wavelength is then 
the wavelength corresponding to the mode that is most 
unstable. In the calculation above, we have picked out 
the mode that in mean field theory corresponds to the 
slowest decaying mode as the wavelength of the quasi- 
patterns. This is because the denominator of the power 
spectrum is equal to the product of the eigenvalues of the 
stability matrix squared. This product is smallest for the 
slowest decaying mode, which is also the mode that will 
go unstable in mean field theory first as parameters are 
varied. Therefore the wavelength of the quasi-patterns 
corresponds to the wavelength of the mean field patterns. 

C. Period of quasi-cycles 

A similar calculation to the calculation above for the 
wavelength of the quasi-patterns can be carried out for 
the period of the quasi-cycles. Consider the denominator 
of the power spectrum with k = 

(pbiP-u: 2 ) 2 +u 2 ((e-p)^) 2 (51) 

Analogous to the wavelength calculation, we seek the 
minimum in u). Simple calculation yields a period of 



Wm \Jlbpi\) — (e — p) 2 ip 2 

Similar arguments to those for the wavelength indicate 
that the period for the quasi-cycles is approximately the 
period for the stable spirals present in mean field theory 

V. DISTINGUISHING QUASI-PATTERNS AND 
QUASI-CYCLES FROM OTHER 
SPATIOTEMPORAL PATTERNS 

To distinguish spatiotemporal patterns generated by 
intrinsic noise from those generated by feedbacks alone 
(i.e. mean field patterns) or by extrinsic noise, it is nec- 
essary to develop theoretical predictions that differ for 
each of these cases. Previous work has focused primarily 
on time series data, focusing on problems such as dis- 
tinguishing quasi-cycles from limit cycles |29 j as well as 
the task of simply determining the amount of extrinsic 
versus intrinsic noise in ecosystems |51) . This work has 
confirmed that both extrinsic noise and intrinsic noise 
are important in real ecosystems for populations such as 
temperate songbirds in Norway, and the beetle species 
Tribolium [STU53] and that quasi-cycles are present in 
real ecological time series data[29]. The work also con- 
firms that the importance of intrinsic noise decreases as 
population density increases, in line with the expectation 
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that the scale of intrinsic noise depends on the scale of 
the population density [53] . 

While separate signatures of quasi-cycles and quasi- 
patterns will be discussed below, one common feature 
that distinguishes quasi-cycles and quasi-patterns from 
their counterparts in mean field theory is that they de- 
pend on the concentration of the population being stud- 
ied. To leading order only the fluctuations have patterns, 
implying that the local populations can be written as 
mean value plus fluctuations scaled by the size of a lo- 
cally well mixed region (see below). Thus the amplitude 
of the patterns relative to the mean population size of 
the fluctuation driven patterns will change as the size of 
a locally well mixed area changes, while the relative am- 
plitude of mean field patterns and limit cycles would not 
change. Such a variation of the size of a locally well mixed 
area could presumably be used to detect quasi-patterns 
and quasi-cycles. 

A. Distinguishing quasi-cycles from limit cycles 

Given a population that has oscillatory abundance in 
time, theory indicates that the oscillations can come from 
either quasi-cycles driven by noise or from population 
density dependent feedbacks alone, perturbed by noise 
(mean field cycles). The key difference mathematically 
is that the power spectrum of limit cycles has a pole at 
its frequency while the power spectrum of quasi-cycles 
does not. In the time domain, this means that the cycles 
driven by intrinsic noise have a short correlation time 
while limit cycles have an infinite correlation time. Since 
poles do not exist in real population data due to stochas- 
ticity, finite size populations, measurement error, etc. 
what this means for real data is that there is a sepa- 
ration of scales between the correlation time of limit cy- 
cles and quasi-cycles. This was first pointed out in detail 
by Pineda-Krch et al. [53] ■ These authors also showed 
that wolverine population cycles are likely quasi-cycles, 
while the celebrated lynx-hare cycles from the Hudson 
Bay company's trapping records are most likely limit cy- 
cles US]. 

Other studies of the role of intrinsic noise have fo- 
cused on intrinsic noise contributions compared to ex- 
trinsic noise contributions as a function of local popula- 
tion size [5TJ[S3]. In frequency space, the best frequencies 
to analyze to distinguish the relative importance of noise 
are high frequencies, corresponding to the short timescale 
fluctuations of the system. To extract predictions for the 
case of intrinsic noise, we look at the large uj asymptotics 
of the power spectrum Eq. [46] at k = 

P(k = 0,u) = w»w m (53) 

where u) m is the modal frequency of the quasi-cycles. 
For cycles driven by extrinsic additive noise, we look at 
the same asymptotics for the power spectrum from the 
heuristic calculation, which, as we noted above, can be 



considered as a calculation for extrinsic noise. In this 
case, the asymptotic form is 

P(fc = 0,aj) = ^ 2+ 4 eV V ), c»c m (54) 

where the variance is independent of population den- 
sity and u> m is the frequency of the quasi-cycles. While in 
this case, both the expressions depend on the square of 
population density, the decay in u> has a power of two for 
intrinsic noise, and of four in the case of extrinsic noise. 
Thus the tails can be easily distinguished in real data. 

B. Distinguishing quasi-patterns from mean field 
patterns 

Similar considerations can be applied to quasi- 
patterns. While further study is needed, the finite 
peak in the power spectrum for quasi-patterns indicates 
that quasi-patterns generically have a shorter correlation 
length than mean field patterns, which have a pole in 
their power spectrum at the wavelength of the pattern. 
Thus the techniques outlined above for distinguishing 
mean field limit cycles from quasi-cycles and applied to 
real populations in [25] translate directly into the space 
domain from the time domain. 

For distinguishing unconserved extrinsic noise and in- 
trinsic noise contributions, the asymptotics for short 
wavelength fluctuations can again be derived for the in- 
trinsic and extrinsic noise cases. For intrinsic noise, we 
have 

P(k,u) = 0) = k~ 2 -, fc>/c m (55) 
v 

where k m is the wave vector of the mode of the power 
spectrum. For extrinsic noise, we have 

P(*,w = 0) = ^- r <fO. (56) 

where the variance is independent of population den- 
sity. Like the quasi-cycle case, the scaling in k differs by 
a power of two between the extrinsic and intrinsic noise 
cases. Contrary to quasi-cycles in the previous section, 
the extrinsic and intrinsic noise lead to different powers 
of population density for large k. This provides a useful 
tool for distinguishing between the effects of unconserved 
extrinsic anoise nd intrinsic noise on the formation of pat- 
terns especially if the density of the populations can be 
varied through comparative study of field data in differ- 
ent ecosystems, or through experiments. These consid- 
erations are quite broad, and should qualitatively apply 
to other systems, such as chemical reaction systems in 
which quasi-patterns or cycles may be present, such as 
the Brusselator model of chemical pattern formation [33] . 

Another possibility beyond those considered here is 
noise manifested through stochasticity in the kinetic pa- 
rameters of the system as is common in ecological models 
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[53] . Systematic study of such a model is well beyond the 
scope of this paper, and is an interesting subject for fu- 
ture research. However, a simple model of the effects of 
weak parameter noise on the linearized dynamics shows 
that the tail of the power spectrum will still be domi- 
nated by the effects of extremely weak extrinsic noise so 
that the k~ A tail of the power spectrum is retained with 
very small fc -8 corrections (See appendix one) indicating 
that weak parameter noise is not a qualitatively impor- 
tant factor in the analysis of quasi-patterns. Absent weak 
extrinsic noise, there is no evidence that weak parame- 
ter noise generates quasi-patterns on its own. Rather, 
it seems to only add corrections to the approach to the 
mean field steady state. 

A potential difficulty with the analysis of power spec- 
trum tails is that preliminary numerical study of quasi- 
patterns suggests that the form of the power spectrum 
used above may only be strictly valid near the onset of 
the patterns due to higher order corrections to the mean 
population. Evidence for this claim is discussed below. 
The implications of this possibility for distinguishing dif- 
ferent kinds of patterns is a subject for future research. 



VI. THERMODYNAMIC LIMIT 

To compare to data, we also must be able to estimate 
the conditions under which the fluctuation driven effects 
described above are important. While the considerations 
that follow are mathematically elementary, they are im- 
portant for the analysis of real data and have not always 
been clearly elucidated in the ecological literature, where 
it has sometimes been assumed that intrinsic noise ef- 
fects are only important if the total population of each 
species is small [54]. In fact, the scale of fluctuations 
is governed instead by the population size in a volume 
(indicated by the parameter V in the calculation above) 
sufficiently small that the time to diffuse out of the vol- 
ume is smaller than the typical time between reactions 
per particle. The confusion arises because when space is 
neglected, the organisms are all confined to such a vol- 
ume, so the scale of fluctuations is determined by the 
total population size [2BJ [SSJ . 

The present calculation shows that there are two sepa- 
rate limits in the construction of reaction-diffusion mod- 
els. One of these limits yields a particular kind of mean 
field theory, and the other, corresponding to what would 
traditionally be called the thermodynamic limit in sta- 
tistical physics, does not yield a mean field theory at all. 
Only in d = do these limits coincide. Recall that the 
theory was constructed by creating a lattice of patches, 
each patch of volume V, and then taking the limit of 
an infinite number of patches, and looking at the con- 
tinuum version of the theory The thermodynamic limit 
corresponds to the limit as the number of patches goes to 
infinity, while the mean field limit corresponds to taking 
the volume of each patch, V, to infinity. 

The parameter V can be estimated by noting that the 



time required on average to diffuse out of a volume V ~ 
L d can be estimated to be t ~ L 2 /D with diffusivity given 
by D. If a characteristic reaction rate for single particles 
(possibly depending in complex way on concentrations) 
is R, then the size of a well mixed volume is constrained 
by 



I?_ 1 
D < R 



(57) 



Rearranging, the size of a well mixed volume is estimated 
by 



V 



d/2 



(58) 



Multiplying the size of a locally well mixed volume by 
the smallest local population density will provide an es- 
timate of how far the system is from the mean field limit. 
However, care should be taken because of amplification 
effects as discussed in the following section. Further in- 
formation about how far the system is from the mean 
field limit can be obtained by examining the properties of 
the power spectrum or correlation functions as discussed 
above. 



VII. VALIDITY OF THE LARGE V EXPANSION 
AND THE SCALE OF QUASI-PATTERNS 

The expansion considered above is only strictly valid 
near the onset of quasi-patterns. While in the absence of 
space, the expansion is valid quite generally, leading to 
excellent agreement between theory and simulation for 
the power spectrum [26] . the spatial structures do not 
seem to be as well captured by the expansion deep in the 
quasi-pattern regime. This is probably due to fluctuation 
corrections to the mean field not studied in the current 
paper. This is suggested by simulated trajectories of the 
reaction-diffusion master equation using the exact algo- 
rithm of Gillespie [SB]. The results of this calculation, 
along with the location on the phase diagram simulated 
are shown in fig. [4] 

The calculation indicates that the patterns deep inside 
the quasi-pattern phase are non-perturbative, due to the 
large variance in populations. We expect that the non- 
perturbative corrections to the mean field solutions arise 
at higher order in the expansion. The analytical theory 
above does not predict the power spectrum of these pat- 
terns, but the calculation of wavelength and period above 
are still approximately valid, since they are obtained by 
finding the least stable modes, which are likely still dom- 
inant, even in the non-perturbative regime. 



VIII. EXPLAINING THE FAILURE OF MEAN 
FIELD THEORY 

From the above calculation, as well as related calcula- 
tions ranging from zero dimensional models of ecosystems 
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FIG. 4: The left hand panel is the phase diagram of the Levin-Segel model with intrinsic fluctuations. The vertical axis is the 
ratio of diffusivities, and the horizontal axis is the dimensionless ratio of generic kinetic parameters. Region I is mean field 
pattern formation, region II is fluctuation driven pattern formation, and region III has no spatial pattern formation. The black 
arrow has its tail on the approximate location in parameter space simulated to produce the spatial patterns shown on the right. 
The right hand panel is a heat map of population density in two dimensions. Note that the number of organisms is highly 
variable, even though mean field predicts no spatial patterns. The fluctuation effects are large, with patch populations ranging 
from 1200 to 0. The axes are the lattice index from simulation. 



to models of biochemical oscillations [2TJ OH] it is 

clear that in many applications where the fundamental 
physics contains intrinsic noise, mean field theory fails 
to describe the oscillatory dynamics in time and space of 
the system even for relatively large systems with many 
degrees of freedom far from a critical point. Qualita- 
tively, this failure can be understood quite generally by 
considering the nature of mean field theory. 

While there are many ways to derive mean field theo- 
ries [55], to understand the failure of mean field theory, 
the simplest approach for systems described by a master 
equation is to note that there are two essential steps to 
deriving a mean field theory: averaging and neglecting 
correlations. 

Consider the first step, averaging. The average of the 



trajectories is given by 

<p = {N(t t x))=Km ^Y, N dt,x) (59) 

where C is the index for realizations of the discrete 
Markov process for the population dynamics, is the 
number of realizations sampled, and N^(t,x) is the real- 
ization of the discrete Markov process. Each individual 
realization may be oscillatory, but the oscillations will 
have a great deal of noise in their amplitude and phase. 
Summing over these oscillatory contributions will under 
many conditions lead to an average, tp, that is no longer 
oscillatory because the variance in amplitude and phase 
between different realizations £ of the stochastic process 
will lead to cancellations of the oscillatory parts in the 
sum for the average above. That is, the sum of noisy 
oscillations is not always oscillatory. Since mean field 
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theory considers the dynamics of averages, it will not 
capture the oscillations present in individual realizations 
of the dynamics unless the feedbacks that generate the 
oscillations are much more important than fluctuations 
(see fig. 






FIG. 5: Sample trajectories of the Markov process for 
predator-prey dynamics. Note that while each is roughly os- 
cillatory, a mean field theory derived from the average of many 
such trajectories would not contain oscillations. 



IX. APPLICATION TO FIELD DATA AND 
EXPERIMENTS 



While the calculation above was intended primarily to 
shed light on the broad theoretical question of the fine 
tuning problem in Turing instabilities rather than the 
Levin-Segel model alone, it would still be satisfying to 
match the predictions above to plankton data. Such an 
application to current field data in planktonic systems 
is very difficult. In part this is because data on plank- 
ton patterns are primarily gathered for large scale spa- 
tial patterns that are driven by turbulent stirring, rather 
than biological interactions as in the theory presented 
here [2]. Convection accounts for most of the spatial 
heterogeneity of plankton at scales above tens of meters 
[T4] . However, there do exist some limited data on plank- 
ton population heterogeneity at meter and shorter length 
scales [13] . Further data on the motility of plankton sug- 
gest that the ratio of diffusivities for predator-prey pairs 
is of order 10 [57] . We calculated above that with generic 
parameters, the criterion 47 yields vj\i> 2.48, while the 



Turing condition yields vj^i > 27.8. Under these condi- 
tions, it is likely that some populations have fluctuation 
driven patterns, if the Turing mechanism is responsible 
for the pattern formation. Current data are not, to our 
knowledge, adequate to go much further. 

There are several additional problems with applying 
the current theory to real planktonic systems, even if 
the data were to be much higher resolution. The first is 
that plankton are enormously diverse, with many species 
interacting with many others, and body sizes and behav- 
iors spanning several orders of magnitude [58] [59]. A 
second problem is that the current theory is so simpli- 
fied that there is no clear connection between many of 
the parameters in the model and what is measured in 
real populations. The best way to carry out the iden- 
tification of quasi-patterns is probably not to engage in 
detailed modeling of the population dynamics, but rather 
to use model independent predictions, such as the den- 
sity dependence of the correlations described above, and 
the power of k and u> for large values of k and u> in the 
power spectrum. Data sets associated with plant systems 
are likely to be amenable to such analysis [IS] . Addition- 
ally, laboratory experiments in engineered microbial [7J, 
or even in chemical systems (see above comments on the 
thermodynamic limit) may provide more controlled ways 
to detect quasi-patterns. 



X. CONCLUSIONS AND PROSPECTS FOR 
FUTURE RESEARCH 

We conclude by noting that our analysis of the model 
in Eqs. [I] has demonstrated that Turing patterns are 
much more generic than is to be expected on the basis 
of mean field theory, partial differential equation anal- 
ysis. We also have pointed out some possible ways in 
which the fluctuation driven spatiotemporal patterns dis- 
cussed can be identified in real data. While this paper 
focused on a single model, we wish to emphasize again 
that the model was deliberately chosen to be generic with 
the goal of providing broad insight into the statistical 
mechanics of the Turing mechanism that can be widely 
applied. As noted in the introduction, the conjectured 
wide applicability of this result has received some sup- 
port from calculations on the Brusselator model [31] and 
a model of Turing patterns in neuronal networks stud- 
ied in the following paper [T7]. Further applications of 
this theory are potentially as wide as the applicability 
of the Turing mechanism, which, as was pointed out in 
the introduction, has been used to explain patterns in an 
enormous variety of systems. In fact, we conjecture that 
perhaps many or most observed Turing patterns are the 
quasi-patterns predicted in this paper. To demonstrate 
this conjecture, the next step is to apply the concepts in 
this paper to an experimentally well-characterized sys- 
tem, such as an engineered bacterial system with Turing 
feedbacks. Another important avenue of investigation 
is to further explore ways to distinguish between quasi- 
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patterns and mean field patterns. Further theoretical 
progress may also be made by addressing with a simi- 
larly detailed theory other noise driven spatiotemporal 
patterns such as intrinsic noise driven epidemic waves, 
which seem to be present in measles and dengue fever 
epidemics [60j 1ST]. 
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Appendix A: Influence of parameter noise on 
quasi-patterns 

To calculate the effects of parameter noise on Turing 
models to a first approximation, we take the linearized 
mean field dynamics 

— iiox = Ax (Al) 

and generalize to 

— iuix = Ax + et;x + Xq (A2) 



where Xo is an initial deviation from mean field equilib- 
rium, e is a small parameter, and £ is a diagonal matrix 
of white noise, variance one. Only contributions to the 
power spectrum that are independent of Xo persist in the 
long time limit. 

Rearranging Eq. |A2| yields 

x = -(A-e£ l y 1 x Q (A3) 

For small e, this can be expanded to yield 

x = A' 1 (1 + eA- 1 ^ x + 0(e 2 ) (A4) 

Taking the dot product and averaging yields the sum of 
the power spectrum for the predator and prey species, 
which is sufficient for detecting the presence or absense 
of quasi-patterns 

(xx*) = (A -1 sc ) (A^xo)* +e 2 A- 2 x {A- 2 x )* (A5) 

Note that this power spectrum depends term by term 
on the initial conditions Xq, indicating that the power 
spectrum is dominated completely by the effects of tran- 
sient patterns at times shorter than the relaxation time 
to the steady state. Parameter noise only has the effect 
of correcting the approach to steady state. Any quali- 
tatively important effects of parameter noise on quasi- 
patterns are thus present only in a nonlinear analysis. 
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